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HIGHLIGHTS 


•  Orbital-free  density  functional  theory  (OFDFT)  is  used  to  study  Li-Si  alloys. 

•  Accurate  cell  lattice  vectors,  bulk  moduli,  and  densities  are  predicted  by  OFDFT. 

•  OFDFT  elastic  constants  and  alloy  formation  energies  agree  well  with  KSDFT. 

•  OFDFT  Li  atom  adsorption  energies  on  a  Si(100)  surface  are  close  to  KSDFT  values. 

•  Linear-scaling  OFDFT  appears  quite  promising  for  large-scale  Li-Si  simulations. 
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Li-Si  interactions  are  of  great  interest  currently  due  to  the  potential  use  of  silicon  anodes  in  Li-ion 
batteries.  As  a  first  step  toward  eventual  nanoscale  characterization  of  lithiation  of  silicon,  here  we 
study  the  crystalline  Li-Si  alloys  LiSi,  Li^Si*  Li7Si3,  Li^U,  Lii5Si4,  and  Li22Sis  using  orbital-free  density 
functional  theory  (OFDFT).  The  recently  proposed  Wang— Govind— Carter  decomposition  (WGCD)  and 
Huang— Carter  (HC)  kinetic  energy  density  functionals  (KEDFs)  are  used  to  evaluate  the  electron  kinetic 
energy.  Both  KEDFs  predict  accurate  cell  lattice  vectors,  equilibrium  volumes,  bulk  moduli,  and  ground- 
state  densities  when  compared  to  Kohn— Sham  density  functional  theory  (KSDFT)  benchmarks.  Elastic 
constants  and  alloy  formation  energies  calculated  with  the  WGCD  KEDF  also  agree  reasonably  well  with 
KSDFT.  Finally,  Li  atom  adsorption  energies  on  the  Si(100)  -2x1  surface  are  calculated  as  a  simple 
initial  test  of  the  Li— Si  mixing  process  during  lithiation  of  silicon.  The  OFDFT  adsorption  energies  again 
are  fairly  close  to  KSDFT  values.  The  results  in  this  work  demonstrate  the  accuracy  of  the  WGCD  and  HC 
KEDFs  for  materials  with  mixed  covalent-metallic  character  and  their  considerable  transferability  under 
different  chemical  environments.  Because  of  its  quasilinear  scaling,  coupled  with  the  level  of  accuracy 
shown  here,  OFDFT  appears  quite  promising  for  large-scale  simulation  of  such  materials  phenomena. 

©  2013  Elsevier  B.V.  All  rights  reserved. 


1.  Introduction 

Lithium  ion  batteries  (LIBs)  feature  many  attractive  advantages 
including  high  energy  density,  light  weight,  and  few  environmental 
concerns,  which  have  led  to  LIBs  dominating  the  rechargeable 
battery  market  in  portable  electronics  in  recent  decades.  At  present, 
most  commercial  LIBs  use  graphite-based  materials  as  the  anode. 
They  are  outstanding  in  terms  of  low  cost,  good  durability,  non¬ 
toxicity,  and  reasonable  energy  density.  However,  to  extend  the 
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use  of  LIBs  in,  e.g.,  storage  for  stationary  power  applications  re¬ 
quires  further  improvement  of  the  energy  density,  which  has 
encouraged  intensive  research  efforts  in  recent  years.  Si-based 
materials  have  been  recognized  as  a  promising  candidate  for  the 
anode  material  [1-4]  due  to  silicon’s  extraordinarily  high  charge 
capacity,  3572  mAh  g~\  compared  to  372  mAh  g  ’  for  graphite  [5- 
7],  However,  during  the  lithiation  process,  the  Si-based  material 
undergoes  significant  structural  change  and  volume  expansion  up 
to  300%  [8—10],  Various  crystalline  phases  including  Lii2Si7,  Li7Si3, 
Li  1 3S14,  LiisSLi,  and  L^Sis  have  been  observed  experimentally 
during  high  temperature  lithiation  [11-17],  At  room  temperature, 
lithiation  also  frequently  induces  a  crystalline-to-amorphous 
transformation  in  the  material  [1—3],  At  full  lithiation  of  silicon, 
the  crystalline  LiisSU  phase  forms  [2,3],  The  complex  structural 
changes  along  with  the  massive  volume  expansion  lead  to  large 


J.  Xia,  EA.  Carter  /  Journal  of  Power  Sources  254  (2014)  62—72 


stresses  in  the  material,  which  cause  particle  fracture,  capacity  loss, 
and  thus  limited  charge/discharge  cycling  capability.  It  is  therefore 
critical  to  understand  the  intrinsic  mechanical  properties  and 
mixing  mechanisms  during  lithiation  and  delithiation  in  Li-Si  al¬ 
loys.  Armed  with  such  insight,  it  may  be  possible  to  render  Si-based 
materials  practical  for  LIB  anode  applications. 

A  number  of  investigations  have  been  carried  out  on  Li— Si 
materials  in  recent  years  [18-28],  In  particular,  first-principles 
computations  using  Kohn-Sham  (KS)  density  functional  theory 
(DFT)  have  provided  many  valuable  predictions  of  mechanical 
properties  and  electronic  structures  for  both  crystalline  and 
amorphous  Li-Si  alloys  [18-21],  Some  KSDFT  studies  [22-28] 
also  investigated  the  dynamics  of  Li  in  Si  anodes  and  mixing 
mechanisms  during  lithiation  or  delithiation.  However,  KSDFT’s 
expensive  (typically  cubic  scaling)  computational  cost  limits 
current  simulations  to  ~200  atoms.  It  thus  prevents  probing 
many  interesting  features  and  phenomena,  which  require  large- 
scale  simulation  (thousands  of  atoms),  such  as  fracture  at  sur¬ 
faces  or  crystalline/amorphous  interfaces  in  nanoparticles  [29— 

31  J. 

Orbital-free  (OF)  DFT  [32],  an  alternative  within  the  DFT 
formalism,  uses  the  electron  density  as  the  sole  variable  and 
computes  energies  using  only  density  functionals.  By  avoiding  k- 
point  sampling  and  orbital  orthogonalization  procedures  necessary 
in  KSDFT,  OFDFT  scales  quasi-linearly  with  system  size  with  a  small 
prefactor.  Calculations  with  thousands  of  atoms  are  now  routine  in 
OFDFT  simulations,  and  a  benchmark  simulation  of  more  than  a 
million  atoms  was  carried  out  more  than  four  years  ago  with 
reasonable  computer  resources  (<200  processors)  [33].  A  number 
of  OFDFT  studies  have  been  reported  in  recent  years,  showing 
instructive  and  promising  results  for  many  practical  materials  and 
properties  [33-44], 

However,  OFDFT’s  enormous  simplification  in  formalism  also 
limits  its  accuracy  and  generality.  Approximating  the  exact  non¬ 
interacting  kinetic  energy  via  kinetic  energy  density  functionals 
(KEDFs)  is  the  major  challenge  in  OFDFT.  Modern  nonlocal  KEDFs 
such  as  the  Chacon-Alvarellos-Tarazona  (CAT)  [45-48],  Wang- 
Teter  (WT)  [49-52],  and  Wang-Govind-Carter  (WGC)  [53]  KEDFs 
all  rely  on  the  Lindhard  linear  response  theory  [54—56]  and 
therefore  are  expected  to  be  accurate  only  for  nearly-free- 
electron-like  systems  such  as  main  group  metals  (since  the  Lind¬ 
hard  response  function  is  the  response  function  of  a  perturbed 
uniform  electron  gas).  Consequently,  the  OFDFT  applications 
mentioned  earlier  are  mostly  confined  to  main  group  metals  and 
their  alloys.  For  other  systems,  such  as  semiconductors,  transition 
metals,  or  molecules,  OFDFT  has  been  far  less  reliable  due,  until 
recently  (vide  infra),  to  the  lack  of  KEDFs  suitable  for  describing 
them. 

The  electronic  structures  of  Li— Si  alloys  are  complicated. 
Experimentally,  many  of  these  alloys  such  as  Lii2Si7  [57],  Li7Si3  [58], 
and  Li15Si4  [59]  have  been  observed  to  be  semiconductors  with 
small  band  gaps.  Nevertheless,  recent  KSDFT  studies  [20,21  ]  pre¬ 
dicted  all  these  alloys  to  be  metallic,  which  is  not  surprising 
considering  the  typical  band  gap  underestimation  by  KSDFT. 
Despite  their  metallic  nature  predicted  by  KSDFT,  the  electron 
density  distributions  in  these  systems  are  still  far  from  nearly-free- 
electron-like.  A  close  examination  of  their  atomic  structures  reveals 
Si— Si  bonds  remaining  in  some  of  the  crystalline  and  also  all 
amorphous  phases  [20,22,27],  where  electrons  are  comparatively 
localized.  Consequently,  these  systems  feature  mixed  metallic- 
covalent  character  and  thus  we  would  expect  that  the  KEDFs 
mentioned  above  will  not  be  able  to  describe  them  properly; 
indeed  this  will  be  shown  in  the  present  work. 

Only  very  recently  have  several  models  been  proposed  aiming  to 
design  advanced  KEDFs  for  covalent  materials  [60,61],  Specifically, 


the  Huang— Carter  (HC)  KEDF  [60]  focuses  on  linear  response 
properties  in  semiconductors,  and  the  Wang— Govind— Carter- 
decomposition  (WGCD)  [61]  model  decomposes  the  total  density 
into  localized  and  delocalized  parts  to  treat  them  separately  with 
different  KEDFs.  They  both  exhibited  remarkable  accuracy  for 
various  semiconductors  [60,61]  and  molecules  [61,62],  However, 
the  materials  examined  in  these  studies  such  as  cubic  diamond 
(CD)  Si,  cubic  zinc  blende  (ZB)  and  hexagonal  wurtzite  III — V 
semiconductors  were  still  relatively  simple  structures.  Further¬ 
more,  both  KEDFs  contain  parameters  that  may  need  adjustment 
depending  on  the  type  of  materials  under  investigation.  Therefore, 
their  transferability  beyond  reported  simple  semiconductors  still 
requires  further  verification  for  more  complex  materials.  In  this 
work,  we  use  the  HC  and  WGCD  KEDFs  to  study  the  crystalline  Li-Si 
alloys  including  LiSi,  Lii2Si7,  Li7Si3,  LfpSU,  LijsSLi,  and  L^Sis.  We 
employ  universal  parameters  in  both  KEDFs  for  a  wide  range  of  Li 
concentrations  in  these  alloys,  aiming  to  test  their  transferability 
under  different  chemical  environments.  OFDFT  using  either  KEDF 
reproduces  reasonably  well  a  variety  of  KSDFT  benchmark  prop¬ 
erties,  including  cell  lattice  vectors,  equilibrium  volumes,  bulk 
moduli,  elastic  constants,  ground-state  densities,  and  surface 
adsorption  energies.  The  results  of  this  work  validate  the  physics 
contained  in  both  KEDF  models,  and  also  promise  a  bright  future  for 
OFDFT  in  modeling  practical  materials  and  properties  beyond 
simple  metals  and  alloys. 

In  the  following,  we  first  briefly  introduce  the  OFDFT  formalism 
and  KEDF  models  in  Section  2.  Numerical  details  are  described 
thereafter  in  Section  3.  In  Section  4,  we  first  discuss  the  use  of  local 
electron-ion  pseudopotentials,  and  then  compare  OFDFT  and 
KSDFT  predictions  of  various  properties  for  different  alloys.  Finally, 
conclusions  and  outlook  are  given  in  Section  5. 

2.  Formalism 

The  Hohenberg— Kohn  theorems  [63]  prove  that  the  electronic 
total  energy  can  be  expressed  as  a  functional  of  the  total  electron 
density.  Specifically,  it  can  be  decomposed  in  the  following  way 
[64]  (all  formulas  use  atomic  units): 

E\fk Otall  =  7's[Ptotal]+/[ftotal]+£xc[Ptotal]+  j  Vext(r)ptotal(r)dr, 

(1) 

where  Ts  is  the  non-interacting  electron  kinetic  energy,  J  is  the 
Hartree  electron— electron  repulsion  energy,  EXc  is  the  electron 
exchange-correlation  (XC)  energy  which  accounts  for  all  errors 
caused  by  approximating  the  exact  interacting  kinetic  energy  and 
electron-electron  interaction  energy  with  Ts  and  J,  Vext  is  the 
external  potential  such  as  ionic  potentials,  and  ptotai  is  the  total 
electron  density.  The  Hartree  and  XC  terms  are  trivial  to  compute  as 
long  as  only  density-dependent  XC  functionals  such  as  local  density 
approximations  or  generalized  gradient  approximations  (GGAs)  are 
used.  However,  the  ionic  potentials  and  Ts  terms  also  need  to  be 
approximated  in  OFDFT.  Evaluating  the  kinetic  energy  with  the 
Laplacian  operator  is  impossible  without  orbitals  and  the 
commonly-used,  more  accurate,  flexible  nonlocal  pseudopotentials 
(NLPSs)  cannot  be  used  in  OFDFT,  also  because  of  the  lack  of  or¬ 
bitals.  Many  local  pseudopotentials  have  been  proposed.  In  this 
work,  we  use  bulk-derived  local  pseudopotentials  (BLPSs)  [62,65— 
67], 

The  major  source  of  error  in  OFDFT  lies  in  the  approximation  of 
Ts  via  KEDFs.  In  this  work,  we  apply  the  WT,  WGC,  CAT  (Sec.  IV  A  in 
Ref.  [48]),  HC,  and  WGCD  KEDFs  to  Li-Si  alloys.  In  modern  KEDF 
theories,  nonlocal  forms  are  usually  adopted  and  the  Ts  term 
generally  can  be  written  as  [52]: 
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Ts  [Ptotal]  =  %[Ptotal] +Tvw[Ptotal] +INL[Ptotal]> 


(2)  Ts[pt otal]  =  (Yf  mil°“‘[ptotal]  -  Tfmil0Cal[Pdel])  +  T™GC[p del], 


where  TTF[ptQtel]  =  3/10(3*)2/3  Jp^dr  [68-70],  TvW[ftotai] 
=  l/8/|VPtotal|2/Ptotaidr  (vW)  [71],  and 


TklPtotai]  =  c/  /  Ptotal(r)w(?(r,r'),  |r-  r'|)p?otal(r')drdr',  (3)  Tf"[p]  =  aTTF[p]  +  bTvW[p], 


(7) 


where  C,  a,  and  /?  are  model-dependent  constants.  The  nonlocal 
term  of  the  CAT  KEDF  [48]  is  a  bit  different,  but  fundamentally 
similar  in  physics.  The  kernel  «  is  determined  by  enforcing  the 
response  function  to  be  the  Lindhard  response  function  when  a 
uniform  electron  density  is  reached.  In  Equation  (3),  the  Fermi 
wave  vector  f(r.r')  can  either  be  a  constant  as  in  the  WT  functional, 
or  density-dependent  as  in  the  CAT  and  WGC  functionals. 

The  HC  KEDF  also  adopts  the  nonlocal  form  shown  above,  but 
additionally  considers  a  gradient  term  to  approximately  account  for 
l/|r  —  r/|  behavior  in  the  response  function  in  semiconductors  (see 
Ref.  [60]  for  the  detailed  formalism).  Specifically, 

f  (r,  i")  =  [37rptotal(r)]1/3  [1+7  ]VPtotf%) .  (4) 

V  Ptotai(r)4/V 

There  are  two  parameters  in  the  model.  A  is  the  parameter  con¬ 
taining  the  key  physics:  a  nonzero  A  takes  the  gradient  term  or  1/ 
|r  —  r'|  into  account  and  therefore  is  optimal  for  semiconductors; 
for  metallic  phases,  A  should  be  set  to  smaller  values  or  zero,  which 
leads  to  a  single-density-dependent  form  of  the  WGC  KEDF.  /3  is  the 
other  adjustable  parameter  (a  4-  /?  =  8/3),  however  it  has  less  in¬ 
fluence  on  most  properties  [60], 

The  WGCD  model  [61]  focuses  on  decomposition  of  the  total 
density  into  a  localized  density  component  such  as  in  covalent 
bond  regions,  and  a  delocalized  density  such  as  in  interstitial  re¬ 
gions  of  all  materials  or  everywhere  in  simple  (non-transition) 
metals: 

Ptotai(r)  =  Pioc(r)  +  Pdei(r),  (5) 

where  pioc  and  pdei  are  the  localized  and  delocalized  densities, 
respectively.  Different  KEDF  models  are  chosen  to  separately 
describe  localized  densities,  delocalized  densities,  and  their  inter¬ 
action  energy.  By  using  the  same  semi-local  KEDFs  for  the  localized 
density  and  interaction  terms,  and  the  WGC  KEDF  for  the  delo¬ 
calized  density,  the  total  non-interacting  kinetic  energy  can  be 
approximated  as: 


Table  1 

Space  groups  and  numbers  of  Li  and  Si  atoms  in  unit  cells  for  all  systems  studied  in 
this  work.  The  kinetic  energy  cutoffs  (in  eV)  used  in  ABINIT  KSDFT-BLPS  and  VASP 
KSDFT— PAW  calculations  and  numbers  of  irreducible  k-points  in  all  KSDFT  calcu- 


Systems  Space  group  #Li/Si  ABINIT  VASP  k-Points 

^cutoff  Ecutoff 


Fd-3m  (227) 

Si  141ja( 88) 

i2Si7  Pnma  (62) 

7Si3  P3212( 153) 

13SU  Pbam  (55) 

15SU  I-43d  (220) 

22Si5  F23  (196) 

lm-3m  (229) 
on  Si  surface  - 


0/2 

8/8 

96/56 

42/18 


88/20 

2/0 

2/24 


900 

900 

600 

900 

900 

900 

600 

900 

900 


400  256 

300  126 

300  27 

300  34 

300  125 

300  10 

300  4 

300  250 

400  28 


as  used  in  previous  work  [61  ].  Another  key  idea  in  the  WGCD  model 
is  to  only  make  use  of  density  information  to  identify  where  elec¬ 
trons  are  localized  and  in  turn  to  use  that  information  to  decom¬ 
pose  the  density.  The  WGCD  scale  function  to  calculate  the 
delocalized  density  has  the  form: 

F(t)  =/(ptotai(r)/p3el-m),  (8) 

and 

PdeiW  =  PtotaiW  x  F(r),  (9) 

where  p[jel  is  the  average  of  the  delocalized  density  and  m  is  a  shift 
parameter.  The  form  of/ is  numerically  constructed  as  in  Ref.  [61  ]. 

There  are  three  parameters  in  the  WGCD  KEDF  besides  the 
original  parameters  in  the  WGC  functional.  The  shift  parameter  m 
in  the  scale  function  is  the  critical  parameter,  which  determines  the 
level  of  density  scaling  in  different  systems,  m  equal  to  zero  is 
optimal  and  transferable  for  various  semiconductors:  larger  m 
values  gradually  revert  the  WGCD  model  back  to  the  WGC  KEDF, 
and  therefore  are  optimal  for  metallic  phases.  The  other  two  pa¬ 
rameters,  a  and  b,  in  the  semi-local  KEDF  have  smaller  effects  on 
most  properties  though  absolute  energies  can  change  slightly  [61  ]. 

3.  Numerical  details 

All  OFDFT  calculations  are  performed  with  our  PROFESS  2.0  code 
[72],  Previously  reported  Li  [67]  and  Si  [62]  BLPSs  constructed  with 
the  Perdew— Burke— Ernzerhof  (PBE)  [73]  GGA  XC  functional  are 
used  in  OFDFT  calculations.  We  compare  our  OFDFT— BLPS  results  to 
KSDFT  with  both  BLPSs  and  the  Projector  Augmented  Wave  (PAW) 
method  [74],  where  the  latter  utilizes  nonlocal  electron-ion  po¬ 
tentials  similar  to  NLPSs.  We  perform  KSDFT— PAW  calculations 
with  the  Vienna  Ab  Initio  Simulation  (VASP)  Package  [75-77],  and 
KSDFT-BLPS  calculations  with  the  ABINIT  code  [78],  The  PBE  [73] 
form  of  GGA  is  used  as  the  XC  functional  in  all  calculations  in  this 
work. 

In  all  calculations,  kinetic  energy  cutoffs  for  the  plane  wave  basis 
are  chosen  to  converge  total  energies  to  within  1  meV  per  atom  for 
each  system.  We  use  1600  eV  for  all  OFDFT  calculations  (equivalent 
to  400  eV  in  ABINIT  in  terms  of  grid  size).  Detailed  kinetic  energy 
cutoffs  and  k-point  sampling  information  for  KSDFT  calculations  are 
listed  in  Table  1  for  all  structures  studied  in  this  work,  along  with 
their  symmetry  space  groups  and  numbers  of  atoms  in  unit  cells. 
Fermi— Dirac  smearing  with  a  smearing  width  equal  to  0.1  eV  is 
used  for  all  structures  except  for  bulk  CD  Si  and  Si(100)  -2x1 
surface  adsorption  where  no  smearing  is  used. 

We  test  the  WT,  CAT,  WGC,  HC,  and  WGCD  KEDFs  in  OFDFT 
calculations.  For  the  WT  KEDF,  a  =  p  =  5/6.  We  use  the  specific  form 
of  the  CAT  functional  following  the  equations  (12),  (23),  and  (24)  in 
Ref.  [48],  where  we  set  a  =  3/5,  /3  =  1,  and  7  =  1.4.  WGC  calculations 
are  very  difficult  to  converge  for  these  Li— Si  alloy  systems  due  to  the 
Taylor  expansion  evaluation  (required  to  retain  quasilinear  scaling) 
[53],  We  have  tuned  large  ranges  of  parameters  of  a,  7,  and  the 
density  expansion  center  p*,  and  find  it  can  only  converge  for  one  of 
the  alloys,  LiisSU,  using  the  universally-derived  density  exponents 
a  =  (5  -  V5)/6,  0  =  5/3  -  a,  7  =  3.6,  and  p*  =  p0  [52,53],  For  the  HC 
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KEDF,  we  adjust  A  and  fl  values  to  obtain  reasonable  bulk  properties 
and  energetics.  We  find  intermediate  values  of  A  =  0.006  and 
=  0.65  produce  optimal  results  for  equilibrium  volumes  and  bulk 
moduli  of  all  alloys  compared  to  KSDFT  benchmarks.  This  param¬ 
eter  set  is  used  to  calculate  all  phases  to  examine  its  transferability. 
In  WGCD  calculations,  we  keep  m  =  0,  which  is  optimal  for  CD  Si,  as 
reported  previously  [61  ].  We  only  adjust  a  and  b  values  in  the  semi¬ 
local  KEDF.  First,  a  =  0.835  and  b  =  0.679  (universal  averaged  values 
from  various  semiconductors  found  in  previous  work)  [61]  are 
tested  to  see  their  transferability  in  these  more  complicated  Li— Si 
alloy  phases.  We  then  slightly  re-adjust  these  two  parameters  and 
find  a  =  0.810  and  b  =  0.6  generate  slightly  better  alloy  formation 
energies  (by  ~0.1  eV)  and  similar  bulk  properties.  These  two 
parameter  sets  are  both  applied  to  calculate  all  phases.  Further  fine 
tuning  can  lead  to  optimal  parameters  for  each  phase,  respectively, 
giving  even  better  bulk  properties  (by  ~  1  %  for  equilibrium  volumes 
and  ~  5  GPa  for  bulk  moduli);  the  optimal  parameters  for  each  alloy 
are  not  identical  but  all  very  close  to  those  above.  Since  the  uni¬ 
versal  parameter  sets  also  predict  reasonable  results  and  using  the 
same  Hamiltonian  for  different  systems  allows  us  to  make  rigorous 
comparisons,  we  do  not  report  the  results  calculated  with  different 
optimal  parameters.  Finally,  the  parameters  in  the  WGC  KEDF  when 
performing  WGCD  calculations  are  the  same  as  reported  previously 
[61],  When  using  the  WGCD  model  to  do  Si  surface  calculations, 
particular  care  needs  to  be  taken  when  calculating  average  den¬ 
sities.  See  details  in  Sec.  Ill  of  Ref.  [61], 

For  each  alloy,  we  compute  equilibrium  volumes,  bulk  moduli, 
and  alloy  formation  energies.  All  structures  are  first  fully  relaxed  in 
KSDFT  calculations.  An  atomic  force  tolerance  of  2  x  10-4  hartree— 
bohr 1  (1  hartree  =  27.2114  eV,  1  bohr  =  0.529  A)  is  used  in  VASP 
and  5  x  10-5  hartree-bohr-1  in  ABINIT.  The  default  stress  tensor 
thresholds  in  both  VASP  and  ABINIT  are  used.  In  OFDFT  calculations, 
we  fix  atomic  internal  coordinates  from  KSDFT-BLPS  calculations 
and  optimize  cells  by  manually  scanning  degrees  of  freedom  in 
lattice  vectors  (stress  tensors  for  the  WGCD  KEDF  are  not  available 
yet  in  PROFESS;  atomic  forces  are  available  but  OFDFT  ion  relaxa¬ 
tion  for  these  complicated  alloy  structures  is  still  not  accurate 
enough  and  needs  further  investigation  that  is  left  for  future  work). 
After  obtaining  equilibrium  lattice  vectors,  we  expand  and 
compress  the  cell  isotropically  by  up  to  2%  with  atomic  internal 
coordinates  fixed.  For  HC  KEDF  calculations  only,  a  ±5%  volume 
change  is  adopted  due  to  numerical  convergence  difficulties  during 
interpolation  (a  large  number  of  interpolation  bins  is  needed  for 
convergence  to  1  meV  per  atom;  a  small  number  of  bins  leads  to 
noisy  curves  on  a  fine  energy  scale).  The  energy  versus  volume 
curves  are  then  fit  to  Murnaghan’s  equation  of  state  [79]  to 
compute  bulk  moduli.  The  alloy  formation  energies  are  calculated 
according  to  the  equation: 

Ef  =  EUxSh_x-xEu-(\-x)Esi,  (10) 

where  £uxsi,_x  is  the  total  energy  per  atom  and  x  is  the  atomic 
concentration  of  Li  in  each  alloy.  £u  is  the  total  energy  per  Li  in  the 
body-centered-cubic  (bcc)  phase  and  Esi  is  the  total  energy  per  Si  in 
the  CD  phase,  which  are  the  ground-state  phases  of  each  element 
under  ambient  conditions. 

We  calculate  Cn,  C22,  C33,  C44,  C55,  C66,  C23,  C13,  C12  for  all  phases, 
and  additionally  C26  for  LiSi.  In  each  elastic  constant  calculation,  we 
apply  a  strain  tensor,  e,  to  the  equilibrium  structure: 


(id 


b) 


Fig.  1.  a)  Side  view  and  b)  top  view  of  the  reconstructed  Si(100)  surface  with  a  p(2  x  1) 
geometry.  Three  different  sites  T,  A,  and  B  are  selected  for  Li  atom  adsorption  calcu¬ 
lations.  The  figure  was  produced  with  the  VESTA  software  [80], 


(Sxx  exy  ezx  \ 

exy  eyy  eyz  I ,  (12) 

ezx  Cyz  Czz  ) 

where  ey  are  strain  components  defined  in  Cartesian  coordinates. 
The  strain  energy  resulting  from  the  deformation  can  be  expressed: 

E  =  E0  +  ^eCeT,  (13) 

where  E o  and  Vo  are  the  total  energy  and  volume  of  the  equilibrium 
structure,  e  =  (.exx,eyy,ezzi2eyz,2exz,2exy),  and  C  is  the  stiffness  tensor. 
For  Cn,  C22,  and  C33,  e  is  set  to  be  (5, 0,0, 0,0,0),  (0,5, 0,0, 0,0),  and 
(0,0, 5, 0,0,0),  respectively,  where  deformation  range  5  is  within  ±1%. 
The  resulting  total  energy  versus  5  data  is  fit  to  the  quadratic 
function 

E  =  E0  +  ^C82V0,  (14) 

to  calculate  the  corresponding  elastic  constant.  Similarly,  for  C44, 
C55,  and  C66,  e  is  set  as  (0,0,0,25,0,0),  (0,0,0,0,25,0),  and  (0,0,0,0,0,25), 
respectively  and  fit  to 


where 


primitive  lattice  vectors  and 


E  =  E0  +  2  C52V0. 

To  calculate  C12,  e  =  (5, 5, 0,0, 0,0)  and  then  we  fit  to 


(15) 
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£  =  £o  +  j(C11  +  C22  +  2C12)52Vo. 

Similarly  e  =  (5, 0,5, 0,0,0)  and  e  =  (0,5, 5, 0,0,0),  and 

(16) 

£  =  £0  +  ^(C11  +  C33  +  2C13)52Vo, 

(17) 

£  =  E0  +  ^  (C22  +  C33  +  2C23)52Vo, 

(18) 

are  used  for  C13  and  C23  calculations,  respectively.  To  calculate  C26, 
e  =  (0,5,0,0,0,25)  and  the  equation 

£  =  £0  +  2  (C22  +  4  C66  +  4C26)52Vq, 

(19) 

are  used  to  do  the  fitting.  Due  to  the  convergence  difficulties  with 
the  HC  KEDF,  we  only  employ  the  WGCD  KEDF  in  OFDFT  to  calculate 
elastic  constants. 

Finally,  Li  atom  adsorption  on  the  reconstructed  Si(100)  -2x1 
surface  is  studied.  We  report  results  of  spin-unpolarized  calcula¬ 
tions;  we  also  did  spin-polarized  calculations  and  found  almost  no 
differences  (<10-3  meV  per  atom).  We  use  a  p(2  x  1)  geometry 
with  12  layers  (24  atoms)  for  the  Si(100)  -2x1  surface.  Two  Li 
atoms  are  then  put  in  selected  sites  above  or  beneath  the  surfaces 
(Fig.  1 )  on  both  sides  of  the  slab  to  avoid  creating  a  net  dipole  within 
our  periodic  boundary  condition  calculation.  10  A  of  vacuum  be¬ 
tween  periodic  slabs  is  used  as  a  buffer  so  as  to  minimize  spurious 
interactions  between  periodic  images.  Geometries  are  fully  relaxed 
in  all  KSDFT  calculations,  while  OFDFT  employs  KSDFT-BLPS 
relaxed  geometries  and  only  optimizes  the  electron  density.  The 
adsorption  energies  are  calculated  as 

fad  »  (f surface  i  Li  -  Surface  -  2ELi)  /%,  (20) 

where  £Surface+iJ  is  the  total  energy  of  the  Si  surface  with  two 
adsorbed  Li  atoms,  FSUrface  is  the  total  energy  of  the  pure  Si  surface, 
and  £Li  is  the  total  energy  for  a  single  Li  atom.  In  the  OFDFT 
calculation  for  single  Li  atom,  the  vW  KEDF  is  used  as  it  is  exact  for 
single-orbital  systems  (here  only  the  2s  electron  of  Li  is  treated 
while  the  Is  electrons  are  subsumed  into  the  pseudopotential). 

4.  Results  and  discussion 

4.3.  Validation  of  bulk-derived  local  pseudopotentials 

Only  local  pseudopotentials  (LPSs)  can  be  used  in  OFDFT  due  to 
the  lack  of  orbitals  (recall  that  NLPSs  contain  orbital-based  pro¬ 
jection  operators).  LPSs  apply  an  identical  potential  to  all  electrons, 
in  contrast  to  NLPSs  or  the  PAW  method  that  apply  different  po¬ 
tentials  to  electrons  of  different  angular  momenta.  In  this  work,  we 


Optimized  b/a  and  c/n  ratios  of  lattice  vectors  in  Li-Si  alloys.  KSDFT— PAW  and 
KSDFT-BLPS  results  are  both  given.  The  OFDFT  results  are  calculated  with  A  =  0.006 
and  (3  =  0.65  in  the  HC  KEDF  and  m  =  0,  a  =  0.835,  and  b  =  0.679  in  the  WGCD  KEDF. 
See  Section  3  for  other  numerical  details. 


LiSi  U7Si3  Li12Si7  Li12Si7  Lil3Si4  U13Si4 

c/a  c/a  b/a  c/a  b/a  c/a 


KSDFT— PAW 
KDFT-BLPS 
OFDFT— WGCD 
OFDFT— HC 
OFDFT— CAT 


0.614  2.370 

0.576  2.393 

0.559  2.431 

0.56  2.48 

0.517 


2.299  1.675 

2.324  1.690 

2.378  1.739 

2.40  1.76 


1.904  0.559 

1.910  0.556 

1.929  0.537 

1.92  0.54 


use  previously  reported  BLPSs  [62,67],  without  including  a  non¬ 
linear  core  correction  [81]  for  either  Li  or  Si.  To  ensure  that  our 
BLPSs  are  sufficiently  accurate,  we  start  with  a  set  of  validation  tests 
within  KSDFT.  We  compare  structures  and  properties  predicted  by 
KSDFT  using  BLPSs  benchmarked  against  the  KSDFT-PAW  method, 
to  identify  errors  due  to  the  BLPSs  alone,  before  moving  on  to 
OFDFT  calculations  which  will  have  additional  errors  due  to  the 
KEDFs. 

The  cell  lattice  vector  ratios,  b/a  and  c/a,  in  relaxed  equilibrium 
structures  are  listed  in  Table  2.  Our  PAW  results  are  all  consistent 
with  those  reported  in  Refs.  [18,20],  The  differences  between 
KSDFT-BLPS  and  KSDFT-PAW  ratios  are  generally  very  small,  of 
the  order  of  1%,  except  for  LiSi  which  exhibits  a  deviation  of  6%.  The 
relaxed  internal  coordinates  (not  reported  here)  with  BLPSs  are  also 
almost  the  same  as  with  PAW,  with  the  average  error  being  less 
than  0.002  in  fractional  coordinates.  The  predicted  equilibrium 
volumes  and  bulk  moduli  for  each  alloy  are  shown  in  Tables  3  and  4, 
and  Fig.  2.  Again,  our  PAW  equilibrium  volumes  are  almost  the  same 
as  found  in  Refs.  [18,20],  while  the  PAW-GGA  results  in  Ref.  [19]  are 
generally  larger  than  all  others.  Our  PAW  bulk  moduli  are  very  close 
to  the  values  in  Ref.  [20],  while  the  bulk  moduli  reported  in  Ref.  [18] 
are  a  bit  smaller  for  Li12Si7  but  similar  for  other  alloys  (there  is  one 
minor  technical  difference  between  the  present  paper  and  the 
other  previous  studies:  in  Refs.  [18,20],  ion  positions  were  relaxed 
for  deformed  structures  when  calculating  bulk  moduli,  while  we 
fixed  ion  positions  as  in  the  equilibrium  phases.  Relaxing  ions  will 
lead  to  slightly  reduced  energies  for  deformed  structures,  and 
consequently  we  expect  our  bulk  moduli  to  slightly  exceed  those 
from  the  other  studies,  but  by  no  more  than  5  GPa.  This  is  similar  to 
the  elastic  constant  results  in  Section  4.3).  Overall,  lattice  constants 
and  equilibrium  volumes  obtained  by  BLPSs  and  PAW  are  very  close 
to  each  other  (Table  3);  LiSi  exhibits  the  largest  deviation  ( ~3.5%), 
consistent  with  the  deviation  in  lattice  vector  ratios.  The  trend  in 
equilibrium  volumes  (increasing  with  Li  concentration)  for  the  al¬ 
loys  is  captured  very  well  (Fig.  2a).  More  significant  differences  are 
seen  in  the  bulk  moduli  (Table  4).  The  KSDFT-BLPS  bulk  moduli  for 
the  alloys  are  generally  somewhat  larger  (by  ~  10  GPa)  than  the 
KSDFT-PAW  moduli,  but  the  trend  of  elastic  softening  as  the  Li 
concentration  increases  is  correctly  captured  (Fig.  2b).  Finally,  we 
calculate  alloy  formation  energies  (Fig.  3).  Our  PAW  results  are 
consistent  with  Fig.  1  in  Ref.  [20]  and  Table  2  in  Ref.  [18],  with 
differences  less  than  0.005  eV  atom-1.  BLPS  errors  compared  to 
PAW  values  are  mostly  very  small,  within  several  tens  of  meV;  only 
for  LiSi  again,  KSDFT-BLPS  calculations  predict  a  bit  larger  forma¬ 
tion  energy  by  ~  80  meV. 

Overall,  the  use  of  BLPSs  does  not  lead  to  large  errors  in  various 
properties  predicted  within  KSDFT,  thereby  justifying  their  use  in 
what  follows.  Next  we  use  BLPSs  in  OFDFT  calculations,  comparing 
against  KSDFT-BLPS  benchmarks,  to  isolate  the  errors  in  various 
KEDFs  while  keeping  the  electron-ion  interaction  the  same  in  both 
sets  of  calculations. 

4.2.  Bulk  properties 

We  first  relax  cell  lattice  vectors.  We  find  it  is  not  a  trivial  task  in 
OFDFT.  The  optimal  fa/a  and  c/a  ratios  obtained  using  the  FIC, 
WGCD,  and  CAT  KEDFs  within  OFDFT  are  listed  in  Table  2  for  the 
alloys  that  contain  degrees  of  freedom  in  their  cell  lattice  vectors. 
We  do  not  report  WGC  KEDF  results  because  it  fails  to  converge  for 
these  alloys.  Although  these  alloys  are  all  predicted  to  be  metallic  in 
KSDFT,  the  density  variation  is  still  very  significant  (see  Section  4.5). 
This  is  caused  by  localized  electrons  existing  within  Si-Si  bonds 
(for  LiSi,  Lii2Si7,  L^Sfy  and  IiuSU)  or  around  Si  atoms  (for  LiisSLj 
and  Li22Si5).  Consequently,  the  WGC  KEDF  becomes  both  physically 
and  numerically  unsound,  and  diverges  when  performing  the 
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Table  3 

Equilibrium  volumes  (Vo,  A3)  per  Si  atom  for  CD  Si  and  each  Li-Si  alloy,  and  V0  per  Li  atom  for  I 
X  =  0.006  and  /}  =  0.65  in  HC  calculations;  m  =  0,a  =  0.835,  and  b  =  0.679  are  used  in  WGCD_A  r 
numerical  details.  The  same  parameter  sets  were  used  in  Table  4. 


CD  Si  LiSi  Li12Si7  Li7Si3 


KSDFT— PAW  20.455 

KSDFT-BLPS  20.444 

HC  21.600 

WGCD_A  20.340 

WGCD_B  20.785 

WGC  23.426 

CAT  19.566 


31.372 

30.278 

30.115 

29.084 

28.970 

29.170 


42.872 

43.143 

42.375 

40.916 

40.049 


49.808 

49.921 

49.457 

49.720 

49.280 

46.723 


Li,  calculated  with  KSDFT— PAW,  KSDFT-BLPS,  and  OFDFT.  We  employ 
i  and  m  =  0,  a  =  0.81,  and  b  =  0.6  in  WGCD_B.  See  Section  3  for  other 


Li  nSi4 


66.176 

66.792 

68.006 

67.503 

67.063 

62.402 


Li,5Si4 


74.563 

75357 

73.791 

76.244 

75.446 

67.676 

69.423 


Li22Si5 


81.251 

83.245 

83.044 

84.071 

83.174 

78.851 


bcc  Li 


20.353 

20.204 

19.803 

19.401 

19.045 

20.208 

20.272 


:h  KSDFT— PAW, 


CD  Si  LiSi  U12Si7  U7Si3  U13Si4  U15Si4  Li22Si5 


WGCD_A 

WGCD_B 

WGC 


39 


Fig.  2.  a)  Equilibrium  volume  per  Si  atom  (A3)  for  CD  Si  and  each  Li-Si  alloy,  and  per  Li 
atom  for  bcc  Li,  and  b)  bulk  modulus  (GPa)  for  each  phase,  calculated  with  KSDFT— 
PAW,  KSDFT-BLPS,  and  OFDFT.  7  =  0.006  and  (3  =  0.65  are  used  in  OFDFT— HC  calcu¬ 
lations;  for  WGCD  results,  we  report  two  data  sets  with  m  =  0,  a  =  0.835,  and 
b  =  0.679,  as  well  as  m  =  0,  a  =  0.81,  and  b  =  0.6. 


required  Taylor  series  expansion  evaluation.  The  WT  and  CAT 
KEDFs,  both  of  which  are  also  based  on  the  perturbed  uniform 
electron  gas,  should  be  inappropriate  for  describing  Li-Si  alloys  for 
the  same  reason.  In  fact,  although  the  WT  KEDF  converges  for  all 
systems,  the  energies  are  unbounded,  i.e.,  energies  always  decrease 
as  volumes  increase,  just  as  found  earlier  for  pure  Si  [82],  The  CAT 
KEDF  performs  better  than  the  WT  KEDF  presumably  because  it  has 
a  density-dependent  kernel.  It  converges  for  every  phase  (only  a 
first  order  expansion  is  needed  to  evaluate  the  CAT  KEDF,  unlike  the 
WGC  KEDF  which  requires  a  second  order  expansion)  and  exhibits 
an  energy  minimum  versus  volume  with  fixed  b/a  and  c/a  ratios. 
However,  no  minimum  can  be  found  when  searching  for  optimal 
ratios  for  any  alloy  except  for  LiSi,  i.e.,  the  total  energies  always 
decrease  with  increasing  or  decreasing  b/a  and  c/a  ratios.  The  CAT 
KEDF  optimized  c/a  ratio  of  LiSi  is  listed  in  Table  2,  but  it  is  not 
accurate  when  compared  to  KSDFT-BLPS  calculations.  Therefore, 
because  of  the  convergence  difficulties  and/or  lack  of  clear  minima, 
all  bulk  properties  shown  hereafter  for  the  WGC  and  CAT  KEDFs  are 
calculated  employing  KSDFT-BLPS  optimized  ratios  (except  for  LiSi 
for  the  CAT  KEDF).  By  contrast,  both  the  WGCD  and  the  HC  KEDFs 
predict  quite  accurate  ratios  for  all  the  alloys,  very  close  to  bench¬ 
mark  KSDFT-BLPS  and  even  to  the  KSDFT— PAW  results.  We  only 
report  HC  ratios  to  the  second  significant  digit  after  decimal  point 
due  to  the  interpolation  convergence  difficulties  mentioned  in 
Section  3. 

We  next  report  equilibrium  volumes  and  bulk  moduli  for  each 
phase  (Tables  3  and  4,  and  also  Fig.  2).  Again,  WGC  results  for  Lil5Si4 
and  CAT  results  for  all  alloys  are  not  satisfactory  with  under¬ 
estimated  equilibrium  volumes  and  overestimated  bulk  moduli.  By 
contrast,  the  HC  KEDF  (that  describes  localized  covalent  bonds 


Fig.  3.  Alloy  formation  energies  (eV)  per  atom  (Ef)  for  each  alloy,  calculated  with 
KSDFT— PAW,  KSDFT-BLPS,  and  OFDFT.  X  =  0.006  and  8  =  0.65  are  used  in  OFDFT— HC 
calculations;  for  WGCD  results,  we  report  two  data  sets  with  m  =  0,  a  =  0.835,  and 
b  =  0.679,  as  well  as  m  =  0,  a  =  0.81,  and  b  =  0.6. 
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Fig.  4.  Total  energy 
energy  minimum  as 


re  used  in  OFDFT-HC  calculations;  for  WGCD  re 


well)  delivers  equilibrium  volumes  and  bulk  moduli  much  closer  to 
KSDFT— BLPS  benchmarks.  Likewise,  the  WGCD  model  with  the 
universal  parameter  set  a  =  0.835  and  b  —  0.679  determined  pre¬ 
viously  [61]  demonstrates  good  transferability  for  these  more 
complicated  structures,  predicting  accurate  equilibrium  volumes 
and  bulk  moduli  for  all  alloys;  another  slightly  adjusted  parameter 
set  a  =  0.81  and  b  =  0.6  also  generates  similar  results.  The  absolute 
energies  (not  shown  here)  reveal  some  errors  with  the  universal 
parameters  but  the  overall  trend  and  relative  energies  are  still 
reasonable  for  the  WGCD  KEDF,  as  shown  when  calculating  alloy 
formation  energies  in  Section  4.4.  Some  energy  versus  volume 
curves  are  plotted  as  examples  in  Fig.  4;  the  shapes  of  OFDFT  curves 
agree  very  well  with  KSDFT  ones  over  a  large  deformation  range. 
The  WGCD  KEDF  with  a  =  0.835  and  b  =  0.679  generally  provides 
the  most  accurate  absolute  total  energies  compared  to  KSDFT  for  all 
alloys  (except  for  LiSi).  Statistically,  80%  of  the  HC  and  WGCD  results 
have  errors  in  equilibrium  volumes  within  3%,  corresponding  to 
less  than  1%  error  in  lattice  constants.  The  largest  error  in  lattice 
constants  is  also  below  2%.  The  predicted  bulk  moduli  are  also  fairly 
accurate  compared  to  KSDFT-BLPS  benchmarks  (errors  within 
5  GPa),  except  the  WGCD  KEDF  overestimates  the  LiSi  bulk  modulus 
by  about  15  GPa.  The  general  trends  of  volume  expansion  and 
elastic  softening  through  the  increasing  Li  concentration  are  clearly 
predicted  by  OFDFT  as  plotted  in  Fig.  2  along  with  KSDFT-BLPS 
results.  Overall,  we  see  very  good  agreement  between  OFDFT  and 
KSDFT  data. 

4.3.  Elastic  constants 

We  also  employ  the  WGCD  KEDF  to  calculate  elastic  constants 
for  every  alloy  as  well  as  for  bcc  Li  and  CD  Si  (Table  5).  Our  PAW 


results  overall  agree  with  Ref.  [18],  with  some  values  a  bit  larger  as 
discussed  in  Section  4.1.  However,  we  also  observe  some  large 
differences,  for  instance,  for  C44  of  LiuSE},  Cn,  and  Cn  of  L^Sis.  We 
found  the  results  reported  in  Ref.  [83]  are  also  very  different  from 
Ref.  [18]  but  close  to  ours.  Here,  we  focus  on  comparing  our  KSDFT 
and  OFDFT  results.  The  HC  KEDF  is  not  used  due  to  its  numerical 


Table  5 

Elastic  constants  (GPa)  calculated  with  KSDFT— PAW,  KSDFT-BLPS,  and  OFDFT— 
WGCD  with  m  =  0,  a  =  0.835,  and  b  =  0.679. 


Method_ C„  C22  C33  C44  C55  C66  C12  C13  C23  C26 


CD  Si  KSDFT— PAW 
KSDFT-BLPS 
OFDFT— WGCD 
LiSi  KSDFT— PAW 
KSDFT-BLPS 
OFDFT— WGCD 
Li,2Si7  KSDFT— PAW 
KSDFT-BLPS 
OFDFT— WGCD 
Li7Si3  KSDFT— PAW 
KSDFT-BLPS 
OFDFT— WGCD 
Li13Si4  KSDFT— PAW 
KSDFT-BLPS 
OFDFT— WGCD 
U15Si4  KSDFT— PAW 
KSDFT-BLPS 
OFDFT— WGCD 
Li22Si5  KSDFT— PAW 
KSDFT-BLPS 
OFDFT— WGCD 
bcc  Li  KSDFT— PAW 
KSDFT-BLPS 
OFDFT— WGCD 


160  160  160  99  99 

153  153  153  104  104 

94  94  94  44  44 

105  105  85  56  56 

130  130  102  60  60 

157  157  148  32  32 

86  114  118  45  33 

113  131  134  49  40 

109  146  137  28  36 

84  84  159  28  28 

114  114  171  27  27 

74  74  150  19  19 

80  94  80  28  25 

94  98  89  32  26 

80  71  68  25  17 

51  51  51  32  32 

63  63  63  35  35 

56  56  56  25  25 

56  56  56  42  42 

52  52  52  42  42 

42  42  42  31  31 

20  20  20  15  15 

16  16  16  16  16 


99  55  55  55 

104  56  56  56 

44  100  100  100 

46  18  37  37  6 

46  22  43  43  3 

40  51  49  49  16 

25  4  11  14 

32  12  19  22 

29  32  33  29 

38  9  5  5 

42  30  1  1 

20  35  24  24 

28  9  5  5 

38  15  11  14 

17  24  25  32 

32  23  23  23 

35  29  29  29 

25  30  30  30 

42  20  20  20 

42  30  30  30 

31  35  35  35 

11  14  14  14 

15  14  14  14 

16  14  14  14 
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issues  discussed  in  Section  3.  For  CD  Si,  OFDFT  unfortunately  pre¬ 
dicts  quite  inaccurate  results  compared  to  KSDFT,  consistent  with 
our  previous  study  [61  ].  As  we  increase  the  Li  concentration,  the 
results  become  increasingly  better.  This  is  presumably  because  we 
have  more  delocalized  electron  density  in  the  system  which  can  be 
properly  described  by  the  WGC  KEDF  buried  in  the  WGCD  model. 
For  Li22Si5,  LiisSU,  or  Lii3Si4,  the  elastic  constants  are  in  fair 
agreement  with  KSDFT-BLPS  results  (L45S4  deformation  curves 
are  shown  in  Fig.  5).  The  values  are  most  accurate  for  bcc  Li  because 
the  WGCD  KEDF  reverts  nearly  to  the  pure  WGC  KEDF,  which 
contains  the  correct  physics  for  a  pure  main  group  metal  such  as  Li. 
The  results  are  worse  for  other  alloys  such  as  LiSi  but  the  magni¬ 
tudes  and  relative  ordering  of  elastic  constants  are  still  acceptable 
for  most  cases.  On  average,  OFDFT  produces  an  error  of  ~  14  GPa  for 
all  alloys  compared  to  KSDFT-BLPS  benchmarks. 

4.4.  Alloy  formation  energies 

Next  we  calculate  alloy  formation  energies  with  the  HC  and 
WGCD  KEDFs.  Universal  parameter  sets  are  used  for  all  phases,  with 
our  aim  to  keep  the  same  Flamiltonian  for  all  systems  and  also  to 
test  the  transferability  of  these  two  KEDF  models.  In  Section  4.2,  we 
discussed  absolute  energy  errors  in  WGCD  results;  similar  errors 
exist  in  HC  results  as  well.  In  previous  work  [60,61],  it  was  shown 
that  although  universal  parameters  lead  to  absolute  energy  errors, 


relative  energies  are  quite  reasonable,  i.e.,  a  common  shift  occurs  in 
all  total  energies.  However,  those  tests  were  made  within  very 
similar  structures,  i.e.,  the  CD  Si  and  ZB  group  III— V  semi¬ 
conductors.  Here,  we  calculate  very  different  phases  ranging  from 
the  purely  metallic  bcc  Li  to  purely  covalent  CD  Si,  and  compare 
their  relative  energies.  We  regard  this  as  the  most  demanding  test 
in  this  work,  because  the  parameter  A  in  the  HC  KEDF  and  m  in  the 
WGCD  model  were  shown  to  have  different  optimal  values  when 
treating  different  materials,  especially  for  covalent  and  metallic 
phases  considered  in  previous  work  [60,61], 

The  resulting  alloy  formation  energies  are  displayed  in  Fig.  3. 
Unfortunately,  we  see  unphysical  behavior  in  the  HC  predictions: 
the  alloy  formation  energies  of  some  alloys  (Li13Si4,  Li15S4,  and 
LbzSis)  are  positive  whereas  KSDFT  predicts  them  all  to  be  negative. 
Even  the  three  negative  formation  energies  are  also  not  satisfactory 
since  they  deviate  significantly  from  the  KSDFT  benchmarks,  except 
for  LiSi.  This  illustrates  the  limited  transferability  of  the  HC  KEDF  in 
different  chemical  environments,  though  other  bulk  properties  are 
rather  reasonable.  The  WGCD  KEDF  generally  performs  better  than 
the  HC  KEDF,  at  least  predicting  all  alloys  to  be  energetically 
favorable  with  negative  formation  energies,  consistent  with  KSDFT. 
For  Li7Si3,  Li13Si4,  Li15Si4,  and  Li22Si5,  OFDFT- WGCD  produces  quite 
accurate  formation  energies,  especially  with  a  =  0.81  and  b  =  0.6. 
However,  for  LinSU  and  particularly  LiSi,  the  alloy  formation 


Fig.  5.  Total  energy  per  unit  cell  versus  6  in  elastic  deformation  of  Cn,  C12,  and  C44  for 
LiI5Si4,  as  calculated  with  KSDFT-BLPS  and  OFDFT- WGCD.  a)  The  OFDFT  curves  are 
shifted  to  have  the  same  equilibrium  energy  (5  =  0)  as  KSDFT  to  assist  comparison;  b) 
curves  without  shifting,  m  =  0,  a  =  0.835,  and  b  =  0.679  are  used  in  OFDFT— WGCD 
calculations. 


Fig.  6.  a)  Maximum  and  b)  minimum  densities  in  each  phase  at  its  equilibrium 
structure,  calculated  with  KSDFT-BLPS  and  OFDFT.  X  =  0.006  and  fl  =  0.65  are  used  in 
OFDFT— HC  calculations;  for  WGCD  results,  we  report  two  data  sets  with  m  =  0, 
a  =  0.835,  and  b  =  0.679,  as  well  as  m  =  0,  a  =  0.81,  and  b  =  0.6. 
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Table  6 

Li  atomic  concentrations  x,  average  density  p0  (in  10  2  bohr~3),  and  /5max/po  for  each 
phase  at  its  equilibrium  structure  calculated  with  KSDFT— BLPS  and  OFDFT— WGCD. 
The  maximum  values  of  the  electron  localization  function  (ELF,^)  from  KSDFT— 
BLPS  calculations,  podeilpo  and  pmax//>odei  from  OFDFT— WGCD  calculations  are  also 
listed,  m  =  0,  a  =  0.835,  and  b  =  0.679  are  used  in  WGCD  calculations. 


CD  Si  LiSi  Li12Si7  Li7Si3  Li,3Si4  Li15Si4  Li22Si5  bcc  Li 


KSDFT— BLPS  p0 
OFDFT— WGCD  p0 
KSDFT— BLPS  pmax/Po 
OFDFT— WGCD  pmax/p0 
ELFmax 

OFDFT/WGCD  podd/po 
OFDFT/WGCD  pmJpodel 


0  0.5  0.632  0.7  0.765  0.789  0.815  1 

2.90  2.44  1.96  1.88  1.61  1.52  1.49  0.73 

2.85  2.55  2.07  1.90  1.60  1.52  1.50  0.78 

2.854  2.930  3.847  3.785  4.283  3.844  3.967  1.108 

2.408  2.712  3.502  3.700  4.330  3.946  3.986  1.135 

0.938  0.928  0.912  0.858  0.861  0.819  0.830  0.581 

0.510  0.536  0.536  0.575  0.598  0.604  0.644  0.960 

4.722  5.064  6.532  6.435  7.243  6.532  6.187  1.182 


energies  are  greatly  overestimated  (too  negative),  thus  giving  an 
incorrect  ordering  of  alloy  formation  stability.  This  demonstrates 
the  transferability  of  the  WGCD  KEDF  still  requires  further 
improvement.  Overall,  the  WGCD  model  correctly  favors  alloy  for¬ 
mation  and  thus  may  prove  adequate  for  simulating  lithiation  of 
silicon,  especially  at  high  Li  concentrations  that  are  of  particular 
practical  interest. 

4.5.  Ground-state  densities 

We  next  compare  OFDFT  and  KSDFT  ground-state  densities.  The 
maximum  and  minimum  densities  in  each  alloy  in  their  equilib¬ 
rium  structures  are  plotted  in  Fig.  6.  The  density  peak  in  CD  Si  is 
underestimated  by  both  HC  and  WGCD.  Encouragingly,  the  un¬ 
derestimation  is  significantly  smaller  for  the  alloys.  In  particular, 
the  maximum  densities  from  WGCD  with  a  =  0.81  and  b  —  0.6  are 
very  close  to  those  from  KSDFT  for  most  of  the  alloys.  Other  OFDFT 
KEDFs  predict  a  bit  smaller  (by  ~  10-3  bohr-3)  maximum  densities. 
Errors  in  the  minimum  densities  of  the  alloys  are  very  small,  on  the 
order  of  10  4  bohr-3.  Qualitatively,  the  WGCD  model  tends  to 
overestimate  the  minimum  density,  indicating  over-delocalization 
presumably  due  to  the  influence  of  the  WGC  KEDF  within  the 
WGCD  model.  By  contrast,  the  HC  KEDF  predicts  slightly  smaller 
minimum  densities  compared  to  KSDFT.  The  general  shapes  of 
density  distributions  are  primarily  determined  by  positions  of  the 
nuclei  and  therefore  are  very  similar  in  OFDFT  and  KSDFT  calcula¬ 
tions.  Overall,  though,  OFDFT’s  maximum  and  minimum  densities 
for  the  alloys  agree  quite  closely  with  KSDFT  densities. 

We  can  also  infer  other  physical  properties  qualitatively  through 
densities.  First,  the  OFDFT— WGCD  average  densities  (Table  6)  are 
very  close  to  KSDFT  ones,  indicating  accurate  equilibrium  volumes. 
Secondly,  the  podeilpo  values  in  WGCD  calculations  (Table  6)  are 
qualitatively  consistent  with  the  fact  that  the  system  tends  to  be 
more  metallic  and  have  more  delocalized  electrons  as  the  Li  con¬ 
centration  grows.  This  ratio  generally  increases  as  the  Li  concen¬ 
tration  increases,  implying  a  larger  portion  of  delocalized  electron 
densities.  In  pure  bcc  Li,  up  to  96%  of  the  electron  density  is 


Li  atom  adsorption  energies  (eV)  on  the  Si(100)  -2x1  surface,  calculated  with 
KSDFT— PAW,  KSDFT— BLPS,  and  OFDFT.  For  the  A  and  B  sites,  the  adsorption  energy 
differences  compared  to  the  T  site  are  also  listed.  A  =  0.006  and  (1  =  0.65  are  used  in 
OFDFT-HC  calculations:  we  employ  m  =  0,  a  =  0.835,  and  b  =  0.679  in  OFDFT- WGCD 


£r_ Ea_ Eb_ Ea-Et  Eb-Et 


KSDFT— PAW  -2.231  -1.790  -1.801  +0.440  +0.430 

KSDFT— BLPS  -1.793  -1.649  -1.662  +0.144  +0.131 

OFDFT— WGCD  -2.345  -2.284  -2.287  +0.060  +0.058 

OFDFT-HC  -2.432  -2.302  -2.184  +0.130  +0.248 


delocalized.  This  trend  demonstrates  the  rationale  behind  the 
WGCD  KEDF  is  physically  reasonable  and  justified.  Thirdly,  we  see 
quite  large  maximum  electron  localization  function  (ELF)  [84,85] 
values  in  the  alloys  (Table  6),  significantly  larger  than  in  bcc  Li. 
This  illustrates  that  localized  electrons  also  exist  in  these  alloys, 
presumably  within  Si— Si  bonds  or  around  Si  atoms.  The  Pmax/Po  and 
Pmax/podei  values  in  WGCD  calculations  (Table  6)  confirm  the  exis¬ 
tence  of  localized  electrons;  the  values  in  alloy  phases  are  very 
large,  even  larger  than  in  CD  Si.  The  average  density  becomes 
smaller  due  to  volume  expansion  with  increasing  Li  concentration, 
but  localized  electrons  still  exist,  producing  higher  pmax/po  ratios. 
Finally,  we  observe  the  trend  of  maximum  ELF  decreasing  with 
increasing  Li  concentration,  which  indicates  that  electron  locali¬ 
zation  is  indeed  weaker  with  more  Li.  In  the  WGCD  formalism, 
Ptotai/podei  is  used  as  the  electron  localization  indicator  due  to  the 
lack  of  orbital  and  ELF  information.  Thus,  pmax/podei  can  be  recog¬ 
nized  as  the  equivalent  counterpart  in  OFDFT  of  maximum  ELF  in 
KSDFT.  However,  at  a  quantitative  level,  we  see  inconsistent  trends 
in  these  two  values.  ELFmax  generally  decreases  as  the  Li  concen¬ 
tration  increases  and  in  CD  Si  it  has  its  largest  value.  This  makes 
sense  since  CD  Si  a  purely  covalent  phase,  pmax/podei  values  from 
OFDFT  exhibit  a  different  trend;  it  first  increases  and  then  decreases 
with  more  Li  in  the  alloy,  and  the  value  in  CD  Si  is  actually  the 
smallest  except  bcc  Li.  The  inconsistency  shows  the  limitation  of 
using  only  densities  as  an  electron  localization  indicator  compared 
to  ELF  in  KSDFT.  Although,  the  WGCD  formalism  is  qualitatively 
correct  in  separating  localized  and  delocalized  densities  as  dis¬ 
cussed  above,  the  specific  method  to  determine  levels  of  electron 
localization  in  different  systems  can  be  further  improved.  However, 
the  results  shown  here  are  still  quite  promising  and  we  believe  a 
more  sophisticated  electron  localization  indicator  and  perhaps 
more  elaborate  decomposition  scheme  will  further  advance  the 
accuracy  and  transferability  of  the  WGCD  KEDF  model. 

4.6.  Li  atom  surface  adsorption 

As  a  final  test,  we  calculate  Li  atom  adsorption  energies  on  the 
reconstructed  Si(100)  surface  (Table  7)  as  first  simple  test  of  Li-Si 
mixing  during  lithiation  of  silicon.  We  employ  a  p(2  x  1)  geometry 
where  each  two  Si  atoms  on  the  surface  form  a  buckled  dimer 
(Fig.  1 ).  We  first  put  Li  atoms  in  the  troughs  (site  T)  between  dimer 
rows.  OFDFT  with  both  HC  and  WGCD  KEDFs  predict  accurate  en¬ 
ergies  compared  to  KSDFT,  especially  compared  to  PAW  results 
through  a  fortuitous  cancellation  of  errors  in  the  BLPS  and  KEDF. 
We  then  push  Li  atoms  deeper  into  the  Si  bulk.  We  choose  two 
different  sites  beneath  the  first  layer  of  the  Si  surface.  KSDFT— PAW 
and  KSDFT— BLPS  give  similar  adsorption  energies  for  these  two 
sites,  while  OFDFT  predicts  more  negative  adsorption  energies  with 
both  KEDFs  (Table  7).  As  Li  atoms  beneath  the  surface  distort  the 
lattice,  the  total  energies  will  increase  and  be  higher  compared  to 
above-the-surface  case  (T  site).  The  adsorption  energy  differences 
for  the  A  and  B  sites  with  respect  to  the  T  site  are  also  calculated 
(Table  7).  We  see  a  considerable  discrepancy  between  KSDFT— PAW 
and  KSDFT— BLPS  calculations  mainly  due  to  different  Ft  values. 
BLPS  energy  differences  are  much  smaller  than  PAW  ones,  meaning 
Li  atoms  may  penetrate  more  easily  into  the  Si  bulk.  The  errors  are 
entirely  caused  by  the  BLPS  approximation  and  illustrate  the 
importance  of  non-local  effects.  Both  KEDFs  within  OFDFT  predict 
positive  energy  differences.  The  WGCD  model  produces  relatively 
smaller  differences,  while  the  HC  KEDF  gives  an  excellent  result  for 
the  A  site  and  an  overestimate  for  the  B  site  compared  to  KSDFT- 
BLPS  results.  This  trend  is  consistent  with  the  alloy  formation  en¬ 
ergy  results  in  Section  4.4.  The  WGCD  KEDF  gives  too  negative 
formation  energies  for  small  Li  concentration  cases;  similarly, 
adsorption  energies  here  are  too  negative,  leading  to 
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Table  8 

Average  computational  wall  time  (seconds)  of  KSDFT-PAW  and  KSDFT-BLPS  calculations  using  32  processors,  and  serial  (one  processor)  OFDFT  calculations  with  different 
KEDFs.  Li22Sis  with  108  atoms  in  the  unit  cell  is  used  for  this  set  of  comparisons. 


OFDFT— WT  OFDFT— CAT  OFDFTGCD  OFDFT— HC  KSDFT-BLPS 

Wall  time  94  124  290  4190  8434 


KSDFT-PAW 

5007 


underestimated  energy  differences  (potentially  leading  to  too  large 
a  driving  force  for  lithiation).  In  contrast,  the  HC  alloy  formation 
energy  for  LiSi  is  very  close  to  the  KSDFT-BLPS  result,  indicating  its 
better  accuracy  for  systems  with  lower  Li  concentrations.  As  a 
result,  the  HC  surface  adsorption  energies  are  quite  accurate, 
especially  for  the  A  site.  Generally,  the  OFDFT  results  are  qualita¬ 
tively  reasonable. 

4.7.  Computational  efficiency 

Finally,  we  compare  computational  cost  for  all  methods.  Because 
our  ultimate  interest  is  to  study  mechanical  deformation  mecha¬ 
nisms  that  will  require  relatively  large  sample  sizes,  numerical  ef¬ 
ficiency  is  also  critical.  In  Table  8,  we  list  the  average  computational 
wall  time  for  the  Li^Sis  electron  density  optimization  (108-atom 
unit  cell)  by  OFDFT  with  different  KEDFs  and  also  by  KSDFT  with 
BLPSs  in  ABIN1T  and  PAW  in  VASP.  Parallelization  is  used  with  32 
CPUs  for  all  KSDFT  calculations.  We  see  the  PAW  method  is 
generally  more  efficient  than  the  BLPS  calculations  due  to  signifi¬ 
cantly  reduced  kinetic  energy  cutoffs  that  can  be  used  in  PAW 
calculations.  All  OFDFT  wall  time  is  reported  with  a  single  core 
calculation  to  ensure  a  reasonable  length  of  time  for  the  OFDFT 
calculation  (our  PROFESS  code  is  also  well  parallelized).  The  WT 
KEDF  features  the  smallest  computational  cost  because  it  has  a 
density-independent  kernel  in  its  non-local  KEDF.  The  CAT  nonlocal 
KEDF  is  a  little  bit  (less  than  2  times)  slower  due  to  its  density- 
dependent  kernel.  The  WGC  wall  time  is  not  reported  since  it 
cannot  converge  this  system;  but  in  a  well  converged  calculation,  it 
is  typically  takes  twice  as  long  as  the  corresponding  WT  calculation. 
The  WGCD  model,  which  involves  several  iterations  of  WGC  cal¬ 
culations  (but  each  iteration  is  much  faster  to  converge),  is 
approximately  three  times  slower  than  WT  calculations.  Finally,  the 
HC  KEDF,  which  suffers  increased  computational  cost  due  to  the 
required  interpolation  step  and  related  convergence  issues,  is 
generally  much  slower  than  other  KEDFs.  Considering  the  different 
numbers  of  CPUs  used  in  the  KSDFT  and  OFDFT  calculations,  all 
OFDFT  calculations  are  orders  of  magnitude  faster  than  KSDFT. 
Moreover,  all  KEDFs  scale  quasi-linearly,  including  the  HC  KEDF, 
and  therefore  they  will  be  increasingly  more  efficient  than  KSDFT 
calculations  for  larger  systems.  From  all  perspectives  including 
accuracy,  transferability,  stability,  convergence  issues,  and  numer¬ 
ical  efficiency,  we  consider  the  WGCD  KEDF  model  to  be  particu¬ 
larly  promising  for  large-scale  simulation  of  Li— Si  alloys. 

5.  Conclusions  and  outlook 

Crystalline  Li— Si  alloys  were  investigated  with  OFDFT.  We  first 
used  KSDFT  calculations  to  verify  the  accuracy  of  BLPSs  compared 
to  PAW  calculations.  Optimized  structures,  bulk  and  elastic  moduli, 
and  also  alloy  formation  energies  generally  show  very  small  errors. 
We  then  employed  BLPSs  in  OFDFT  calculations  to  test  various  KEDF 
models.  In  particular,  the  HC  and  WGCD  KEDFs  recently  proposed 
for  covalent  materials  were  closely  examined.  These  two  func¬ 
tionals  predict  quite  reasonable  cell  lattice  vectors,  equilibrium 
volumes,  and  bulk  moduli  for  all  alloys.  They  represent  a  striking 
improvement  over  other  modern  KEDFs  including  the  WT,  WGC, 
and  CAT  functionals,  which  either  fail  to  converge  or  produce 


unphysical  results.  Elastic  constants  computed  with  the  WGCD 
KEDF  are  fairly  reasonable  especially  for  those  alloys  with  higher  Li 
concentrations.  Alloy  formation  energies  were  also  calculated  using 
universal  parameter  sets  in  the  HC  and  WGCD  KEDFs.  The  HC  KEDF 
exhibited  unsatisfactory  transferability  among  most  alloys, 
although  an  accurate  result  was  predicted  for  the  low  Li  concen¬ 
tration  alloy  LiSi.  By  contrast,  the  WGCD  functional  favors  forma¬ 
tion  of  all  alloys  examined  (as  does  KSDFT)  and  generates  good 
formation  energies  for  high  Li  concentration  alloys,  but  greatly 
overestimates  the  formation  energies  of  LiSi  and  Li12Si7.  A  similar 
problem  is  found  when  calculating  Li  adsorption  energies  on 
Si(100),  which  can  be  considered  as  an  extremely-low-Li- 
concentration  case.  The  WGCD  KEDF  underestimates  the  energy 
differences  between  T  and  A/B  sites  while  the  HC  KEDF  gives  larger 
and  more  accurate  results  compared  to  KSDFT  benchmarks.  How¬ 
ever,  both  functionals  give  qualitatively  correct  trends.  Finally,  we 
compared  ground-state  densities  in  each  alloy.  OFDFT  densities  are 
very  close  to  KSDFT  benchmarks,  as  judged  by  comparing 
maximum  and  minimum  densities.  A  closer  examination  of  density 
information  shows  that  both  electron  localization  and  the  fraction 
of  electrons  that  are  delocalized  are  reasonably  captured  within  the 
WGCD  formalism.  However,  the  level  of  electron  localization  is  not 
precisely  determined;  more  sophisticated  electron  localization  in¬ 
dicators  or  WGCD  scale  functions  may  help  to  improve  the  KEDF 
model. 

Overall,  the  OFDFT  results  are  quite  encouraging  with  both  of 
these  nonlocal  KEDFs  designed  for  semiconductors.  With  the 
outstanding  computational  efficiency  shown  in  Section  4.7,  OFDFT 
has  the  ability  to  study  much  larger  sample  sizes  than  KSDFT  is 
capable  of  studying.  It  can  thus  be  used  to  study  thicker  Si  nano¬ 
wires  or  large  surface  structures  while  mimicking  experimental 
conditions  more  realistically.  Bulk  properties,  surface  energies  [61  ], 
and  elastic  softening  are  correctly  reproduced  by  OFDFT,  so  it  may 
be  used  to  investigate  mechanical  deformation  and  failure  during 
lithiation  and  delithiation.  However,  both  KEDFs  are  not  accurate 
enough  to  reproduce  all  alloy  formation  energies,  which  are  closely 
related  to  the  open  circuit  voltage  [86,87],  One  needs  to  be  careful 
when  calculating  related  quantities  within  the  present  OFDFT  ap¬ 
proximations,  especially  with  the  HC  KEDF  which  predicts  many 
positive  alloy  formation  energies.  However,  the  HC  KEDF  is  quite 
accurate  for  low  Li  concentrations,  which  means  that  it  can  be  used 
to  study  the  beginning  stages  of  lithiation.  The  reasonable  HC  Li- 
atom  surface  adsorption  energies  also  indicate  it  would  be  prom¬ 
ising  to  study  Li  atom  diffusion  in  Si  nanowires,  as  studied  in 
Ref.  [25],  By  contrast,  the  WGCD  KEDF  overestimates  alloy  forma¬ 
tion  energies  for  low  Li  concentration  alloys.  The  resulting  larger 
driving  force  of  initial  Li— Si  mixing  might  also  overestimate  the  rate 
of  mixing  during  the  initial  lithiation  process.  However,  it  is  very 
accurate  for  high  Li  concentration  alloys,  starting  from  LiySb,  aside 
from  LiSi.  Consequently,  we  would  use  the  WGCD  KEDF  to  study  the 
latter  part  of  lithiation  process  or  delithiation  processes,  whenever 
the  Li  concentration  is  relatively  high.  Both  KEDFs  are  clearly  not 
perfect  yet.  They  are  not  transferable  enough  for  all  phases,  so 
unfortunately  the  entire  lithiation  and  delithiation  process  cannot 
be  modeled  accurately  with  one  universal  KEDF  yet.  This  is  the 
main  limitation  of  OFDFT  at  present,  and  care  must  be  taken  when 
choosing  a  suitable  KEDF  model  to  explore  the  problem  of  interest. 
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Future  studies  will  work  on  KEDF  models  to  improve  especially 
alloy  formation  energies.  Meanwhile,  we  will  test  OFDFT  accuracy 
on  amorphous  Li-Si  alloys,  since  anode  materials  are  observed  to 
undergo  crystalline-to-amorphous  transformations  at  room  tem¬ 
perature.  Molecular  dynamics  (MD)  simulation  using  OFDFT  forces 
will  be  necessary  as  well.  A  number  of  OFDFT  MD  studies  [67,88— 
91  ]  have  been  reported  in  literature  and  showed  excellent  results. 
Large-scale  OFDFT  MD  simulations  should  be  a  useful  tool  to  probe 
various  interesting  and  critical  questions  in  the  silicon  lithiation 
process,  potentially  offering  new  insights  for  anode  material 
development. 
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